package.list <- c("here", "tidyverse",
"SIBER", "glmmTMB",
"MuMIn", "emmeans",
"ggeffects", "DHARMa",
"effects")
## Installing them if they aren't already on the computer
new.packages <- package.list[!(package.list %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
## And loading them
for(i in package.list){library(i, character.only = T)}
## here() starts at /Users/Ana/Dropbox/2020/Collaborations/dna_isotopes
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## This is DHARMa 0.3.2.0. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa') Note: Syntax of plotResiduals has changed in 0.3.0, see ?plotResiduals for details
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
Load data, which includes cane spider isotope data, plant isotope data (for correcting), cane spider body size data, and island attributes
spider_iso <- read.csv(here("data", "isotopes", "2015 Palmyra Cane Spider Isotopes.csv"))
plant_iso <- read.csv(here("data", "isotopes", "2009-2012 Palmyra Plant Isotope Data.csv"))
spider_size <- read.csv(here("data", "size", "Spider sizes 2009-2015.csv"))
islands <- read.csv(here('data', "Island_data.csv"))
Correct raw \(\delta^{15}N\) values with plant \(\delta^{15}N\)
Calculate SIBER objects per island, including:
\(\delta^{15}N\) range (NR): larger NR, larger trophic diversity within the population
\(\delta^{13}C\) range (CR): larger CR, more resource pools, so more niche diversification among individuals
Total area (TA): niche space occupied by population
Mean distance to centroid (CD): a measure of how similar individuals are to each other that is more invariant to outliers than TA
Nearest neighbor distance (NND): another representation of how similar individuals are to each other in trophic space, or how similar their niches are
Standard deviation in nearest neighbor distance (SDNND): less influenced by sample size than NND, gives a metric of how similar individual niches are to each other
Pre-statistics data visualizations
Statistical analyses (resource use could be any or all measures above of trophic diversity and niche space or raw values for nitrogen or carbon isotopes):
str(spider_iso)
## 'data.frame': 76 obs. of 7 variables:
## $ ID : chr "Castor 1" "Castor 10" "Castor 2" "Castor 3" ...
## $ Amount_ug: num 580 527 439 562 593 ...
## $ d15N : num 19 14.6 17.1 16.8 16.8 ...
## $ Wt...N : num 12.2 12 10.9 11.4 12 ...
## $ d13C : num -23.4 -23.7 -24.1 -24 -22.9 ...
## $ Wt...C : num 42.2 41.5 38.5 39.7 43 ...
## $ Year : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
Need to split ID column into island and ID columns
spider_iso <- spider_iso %>%
separate(ID,
into = c("Island", "ID"),
sep = " "
)
Make sure that island names in spider_iso and plant_iso match
unique(spider_iso$Island)
## [1] "Castor" "Dudley" "Fern" "Leslie" "Lost" "Para" "Sand"
unique(plant_iso$Island.name)
## [1] "Frigate" "Eastern" "Lost" "Paradise" "Holei"
## [6] "Castor" "Fern" "Kaula" "Leslie" "Sand"
## [11] "Engineer" "Portsmouth" "Sacia" "" "Ainsley"
## [16] "Aviation" "Dudley" "N. Fighter" "NS Causeway" "Pollucks"
## [21] "S. Fighter" "Whipoorwill" "Home" "Bird" "Claire"
## [26] "Papala"
‘Para’ needs to be changed to ‘Paradise’ in spider_iso. I’ll do this with an ifelse statement in the mutate function in dplyr, which says “Remake this variable Island such that if Island =”Para“, change it to”Paradise“, otherwise keep the current value in the column for Island”.
spider_iso <- spider_iso %>%
mutate(Island = ifelse(Island == "Para", "Paradise", Island))
There are multiple ways to do this next step, but what I’m going to do is 1) average all the plant \(\delta^{15}N\) values per island giving it the name plant_d15N, and then 2) join that dataframe to the spider dataframe. Last, 3) I can just create a new column with the mutate() function that is the corrected \(\delta^{15}N\) of animal \(\delta^{15}N\) - plant \(\delta^{15}N\).
plant_iso <- plant_iso %>%
group_by(Island.name) %>%
summarise(plant_d15N = mean(d15N))
## `summarise()` ungrouping output (override with `.groups` argument)
spider_iso <- spider_iso %>%
dplyr::select(Island, ID, d15N, d13C) %>%
left_join(plant_iso, by = c("Island" = "Island.name"))
spider_iso <- spider_iso %>%
mutate(d15N_c = d15N - plant_d15N)
(from Layman et al. 2007). the laymanMetrics() function in SIBER calculates the following metrics (explained in ecological context in the roadmap):
\(\delta^{15}N\) range (NR): trophic diversity (higher = greater diversity)
\(\delta^{13}C\) range (CR): resource pool use (higher = broader resource use)
Total area (TA): niche of community (bigger = bigger niche)
Mean distance to centroid (CD): similar to TA, but less influenced by outliers (bigger = bigger niche)
Nearest neighbor distance (NND): trophic similarity (bigger = more individual variation)
Standard deviation in nearest neighbor distance (SDNND): similar to NND but less influenced by sample size (bigger = more individual variation)
Annoyingly, the standard SIBER functions are automatically programmed to take data in a particular format (different species within communities), which is not what we are trying to accomplish here. We are trying to understand relationships between individuals in a population. So - I’m going to use the laymanMetrics() function to manually calculate all of these values for each community (“Island”). Then we can use these in analyses. The laymanMetrics() function asks for two vectors of values which correspond to the two isotope values per individual. What I’m going to do is 1) split the dataframe for isotopes into a dataframe per island, 2) build a function that will compute the layman metrics for each of these split dataframes, and 3) combine those all back into one dataframe so we can compare them.
island_isos <- split(spider_iso, spider_iso$Island)
population_metrics <- function(df) {
x <- df$d13C
y <- df$d15N_c
layman <- laymanMetrics(x, y)
metrics <- as.data.frame(layman$metrics)
metrics <- metrics %>%
rownames_to_column(var= "metric") %>%
rename("value" = "layman$metrics")
return(metrics)
}
island_layman_metrics <- lapply(island_isos, FUN = population_metrics)
island_metrics <- dplyr::bind_rows(island_layman_metrics, .id = "island_layman_metrics")
Let’s attach our island level metadata on size and productivity to this so we will be ready to do analyses! First we’ll rename this column in the metrics dataframe and make sure the naming conventions match between this dataframe and that of the enviornmental data.
island_metrics <- island_metrics %>%
rename("Island" = "island_layman_metrics")
unique(island_metrics$Island)
## [1] "Castor" "Dudley" "Fern" "Leslie" "Lost" "Paradise" "Sand"
unique(islands$Island)
## [1] "Ainsley" "Aviation" "Castor" "Dudley"
## [5] "Eastern" "Engineer" "Fern" "Frigate"
## [9] "Holei" "Home" "Home NE" "Home SW"
## [13] "Kaula" "Leslie" "Lost" "North Fighter"
## [17] "NS Causeway" "Papala" "Paradise" "Pollucks"
## [21] "Portsmouth" "Sacia" "Sand" "South Fighter"
## [25] "Strawn" "Whipoorwill"
Looks good! So we can combine these. And add some categorical values for productivity and size while we’re at it. The ifelse() statements below are saying, look at the values for Island_prod and Island_Area, if this value is above a certain threshold, this is a “high productivity” and/or “big” island, otherwise it is a “low productivity” or “small” island. I’m also going to mutate the order of the metric variable so we can look at the graph below in the order we’ve been talking about these variables.
island_metrics <- island_metrics %>%
left_join(islands, by = "Island") %>%
mutate(prod_level = ifelse(Island_prod > 0.007020000, "high", "low"),
size_level = ifelse(Island_Area > 12339.021, "big", "small")) %>%
mutate(metric = ifelse(metric == "dX_range", "dC_range",
ifelse(metric == "dY_range", "dN_range",
metric))) %>%
mutate(metric = factor(metric, levels = c("dC_range", "dN_range", "TA", "CD",
"NND", "SDNND")))
These data aren’t “tidy” for analyses right now, but they are set up nicely for a quick visualization based on island environmental variables.
ggplot(island_metrics, aes(x = prod_level, y = value, fill = prod_level)) +
geom_boxplot() +
facet_wrap(~metric, scales = "free") +
theme_bw()
### Island size metrics
ggplot(island_metrics, aes(x = size_level, y = value, fill = size_level)) +
geom_boxplot() +
facet_wrap(~metric, scales = "free") +
theme_bw()
Looks like there are stronger patterns by productivity (to be expected based on all the previous literature from this field site), and maybe some patterns with island size, but less clear. We can run some statistics on these down the road, so stay tuned!
#Pre-statistics data visualizations
Questions 1. & 2. above are best asked with the same model (e.g. including both factors and seeing whether they’re important). However, this question encompasses many different values, including:
spider_iso <- spider_iso %>%
left_join(islands, by = "Island") %>%
mutate(prod_level = ifelse(Island_prod > 0.007020000, "high", "low"),
size_level = ifelse(Island_Area > 12339.021, "big", "small"))
ggplot(spider_iso, aes(x = d13C, y = d15N_c, level = prod_level)) +
geom_point(aes(shape = prod_level)) +
stat_ellipse() +
stat_ellipse(aes(color = Island), linetype = "dashed") +
theme_bw()
Values higher on the y in this graph are higher trophic levels (e.g. eating predators vs. herbivores), values farther to the right on the x-axis correspond to more marine resources, more to the right are more terrestrial. It definitely looks like there are some island productivity patterns in bi-plot space!
These same patterns in isotope bi-plot space don’t seem to be happening with islet size.
ggplot(spider_iso, aes(x = d13C, y = d15N_c, level = size_level)) +
geom_point(aes(shape = size_level)) +
stat_ellipse() +
stat_ellipse(aes(color = Island), linetype = "dashed") +
theme_bw()
ggplot(spider_iso, aes(x = prod_level, y = d15N_c, fill = size_level)) +
geom_boxplot() +
theme_bw()
Not at all surprisingly, on high productivity islands, spiders are eating at higher trophic levels, which is what we knew from Young et al. 2013.
ggplot(spider_iso, aes(x = prod_level, y = d13C, fill = size_level)) +
geom_boxplot() +
theme_bw()
In this graph, a value higher on the y corresponds to more marine. Interestingly (to me), spiders on low-productivity islands seem to be taking advantage of more marine resources.
These are the same graphs as above, but as a refresher here.
Island productivity:
ggplot(island_metrics, aes(x = prod_level, y = value, fill = prod_level)) +
geom_boxplot() +
facet_wrap(~metric, scales = "free") +
theme_bw()
Island size:
ggplot(island_metrics, aes(x = size_level, y = value, fill = size_level)) +
geom_boxplot() +
facet_wrap(~metric, scales = "free") +
theme_bw()
Question 3. above would be best addressed with thinking about how body size is related to d15N, d13C or both combined.
do larger spiders have higher d15N?, do larger spiders feed on different resource pools (different d13C)?
Excitingly for this dataset, there doesn’t seem to be a relationship between island size or productivity with spider size. So any answers we get won’t be confused by just island-level differences, which is cool. Also interesting to see these isotope differences among a population that is really similar in body size (IMO)!
spider_size <- spider_size %>%
filter(Year == 2015) %>%
dplyr::select(ID, Length_mm, Mass_g) %>%
separate(ID,
into = c("Island", "ID"),
sep = " "
)
## Warning: Expected 2 pieces. Additional pieces discarded in 41 rows [98, 99, 100,
## 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
## 117, ...].
spider_iso <- spider_iso %>%
left_join(spider_size, by = c("Island", "ID"))
ggplot(spider_iso, aes(x = prod_level, y = Mass_g, fill = size_level)) +
geom_boxplot() +
theme_bw()
Also there don’t seem to be really strong relationships between mass and isotope values
ggplot(spider_iso, aes(x = Mass_g, y = d15N_c)) +
geom_point() +
geom_smooth(method = "lm", se =F) +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
ggplot(spider_iso, aes(x = Mass_g, y = d13C)) +
geom_point() +
geom_smooth(method = "lm", se =F) +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
Or in bi-plot space
ggplot(spider_iso, aes(x = d13C, y = d15N_c, size = Mass_g)) +
geom_point() +
theme_bw()
The differences across islands don’t seem to be driven by body size differences, which I think is super interesting! Instead, they are similar populations using different environments differently!
The statistical analyses we’ll be using for these datasets vary based on the data structure for each data type.
raw isotope values and body size: these data are at the individual level within populations on islands that have similar environmental parameters. Because individuals are grouped in island-level populations (and thus, these individuals are likely to be more similar to each other than other individuals), but what we care about is actually the environmental parameters that describe that group of islands (e.g. how big they are, their base resource productivity), we will use linear mixed effects models (LMM) for these analyses. These models can handle the fact that individuals are grouped by island (a “random effect” in the mixed effect model) while also asking questions about larger variables (a “fixed effect” in the mixed effect model, including island size and productivity values).
layman metrics: these are actually a much easier dataset to deal with, and will only involve using linear models (LM) or generalized linear models (GLM) like “fancy ANOVAs” to ask whether these metrics vary by island productivity and size.
I performed all of these for you and help you with interpretations, however, this is a multi-step and complicated process. I encourage you to download and read this tutorial from my GitHub as you work through the models so you understand why I performed these analyses this way. I’ll just say - these data played very nice - I enjoyed working through these models because they didn’t have many of the issues I usually have with data distributions that make me want to pull my hair out.
To refresh on your questions, which I’m going to go ahead and break down more specifically in the structure of the models I just described
raw \(\delta^{13}C\) ~ productivity*size + (1|Island) (LMM)
raw \(\delta^{15}N\) ~ productivity*size + (1|Island) (LMM)
dC_range ~ productivity*size (LM)
dN_range ~ productivity*size (LM)
TA ~ productivity*size (LM)
CD ~ productivity*size (LM)
NND ~ productivity*size (LM)
SDNND ~ productivity*size (LM)
size ~ productivity*size + (1|Island) (LMM)
raw \(\delta^{15}N\) ~ size + (1|Island) (LMM)
raw \(\delta^{13}C\) ~ size + (1|Island) (LMM)
raw_c <- glmmTMB(d13C ~ prod_level*size_level + (1|Island),
data = spider_iso)
dredge(raw_c)
## Fixed terms are "cond((Int))" and "disp((Int))"
raw_c_best <- glmmTMB(d13C ~ prod_level + (1|Island),
data = spider_iso)
rawC_sim <- simulateResiduals(fittedModel = raw_c_best, plot=T)
summary(raw_c_best)
## Family: gaussian ( identity )
## Formula: d13C ~ prod_level + (1 | Island)
## Data: spider_iso
##
## AIC BIC logLik deviance df.resid
## 185.2 194.5 -88.6 177.2 72
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Island (Intercept) 0.267 0.5167
## Residual 0.506 0.7113
## Number of obs: 76, groups: Island, 7
##
## Dispersion estimate for gaussian family (sigma^2): 0.506
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -24.8169 0.2799 -88.65 < 2e-16 ***
## prod_levellow 1.7160 0.4281 4.01 6.11e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(raw_c_best))
#### raw \(\delta^{13}C\) Interpretation
Predators on high productivity islands have lower \(\delta^{13}C\) values than those on low productivity islands. Lower (more negative) \(\delta^{13}C\) values correspond to more terrestrial resource pools, suggesting that these predators rely more on terrestrial than marine resources. (\(\beta\) = 1.72, p-value < 0.001).
raw_n <- glmmTMB(d15N_c ~ prod_level*size_level + (1|Island),
data = spider_iso)
dredge(raw_n)
## Fixed terms are "cond((Int))" and "disp((Int))"
rawN_sim <- simulateResiduals(fittedModel = raw_n, plot=T)
summary(raw_n)
## Family: gaussian ( identity )
## Formula: d15N_c ~ prod_level * size_level + (1 | Island)
## Data: spider_iso
##
## AIC BIC logLik deviance df.resid
## 245.1 259.0 -116.5 233.1 70
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Island (Intercept) 0.1576 0.3969
## Residual 1.1562 1.0753
## Number of obs: 76, groups: Island, 7
##
## Dispersion estimate for gaussian family (sigma^2): 1.16
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.6001 0.5125 18.731 < 2e-16 ***
## prod_levellow -7.0969 0.7320 -9.695 < 2e-16 ***
## size_levelsmall -2.1764 0.5921 -3.676 0.000237 ***
## prod_levellow:size_levelsmall 2.6411 0.8692 3.038 0.002378 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(raw_n))
Predator trophic position is influenced by a combination of island size and productivity, with predators on high productivity islands having higher \(\delta^{15}N\) values than those on low productivity islands. However, this relationship is mediated by island size, with predators on large high productivity islands having a higher \(\delta^{15}N\) value than those on small high productivity islands. (need to check on beta values and p-values for this model)
First, make the island metrics dataframe wider, meaning there is a column for each metric and its value
island_metrics <- island_metrics %>%
pivot_wider(names_from = "metric",
values_from = "value")
C_range <- lm(dC_range ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(C_range)
## Fixed term is "(Intercept)"
Due to the small sample size, I won’t be surprised if these metrics have lower significance values overall. Based on the dredge() call above it looks like the best model includes no terms. However, because the \(\Delta\)AIC value is close to the cutoff of 2, I’m going to build a model with that variable (prod_level). Likely any significance, if any, will be marginal, in this final model.
C_range_best <- lm(dC_range ~ prod_level,
data = island_metrics)
sim_c_range <- simulateResiduals(fittedModel = C_range_best, plot = T)
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
summary(C_range_best)
##
## Call:
## lm(formula = dC_range ~ prod_level, data = island_metrics)
##
## Residuals:
## 1 2 3 4 5 6 7
## -0.7133 -0.5375 0.5867 0.4125 -0.0875 0.1267 0.2125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7775 0.2628 6.765 0.00107 **
## prod_levellow 0.8658 0.4014 2.157 0.08349 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5255 on 5 degrees of freedom
## Multiple R-squared: 0.482, Adjusted R-squared: 0.3784
## F-statistic: 4.653 on 1 and 5 DF, p-value: 0.08349
Although predators across high-low productivity islands vary in their raw \(\delta^{13}C\) values, they do not significantly vary in the range of values for carbon across the population. This suggests that although across islands they are using different resource pools, within an island population, this does not translate to a wider range of resources used or available to predators. [this could be an interesting discussion point thinking about whether these resources are or are not even available in these different habitat patches or whether predators are just choosing something based on optimal foraging theory]. There was a marginal effect of island productivity on carbon value range (\(\beta\) = 0.87, p-value = 0.08)
N_range <- lm(dN_range ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(N_range)
## Fixed term is "(Intercept)"
That cutoff of \(\Delta\)AIC is higher for this one, so going to go ahead and say that N range does not change with either island variable, even marginally.
N_range_best <- lm(dN_range ~ 1,
data = island_metrics)
sim_n_range <- simulateResiduals(fittedModel = N_range_best, plot = T)
Although predators across high-low productivity and large-small islands vary in their raw \(\delta^{15}N\) values, they do not significantly vary in the range of values for nitrogen across the population. This suggests that they have different trophic niches across different islands, not that island environmental context encourages a broader diversity of trophic resources. This could, again, either be because these resources don’t exist, or some degree of “choice” by predators based on optimal foraging.
TA <- lm(TA ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(TA)
## Fixed term is "(Intercept)"
Again, this \(\Delta\)AIC value is high, so yeah, no marginal patterns to explore.
TA_best <- lm(TA ~ 1,
data = island_metrics)
sim_ta <- simulateResiduals(fittedModel = TA_best, plot = T)
Predators in different environmental contexts do not have broader trophic niches, even though their trophic niche is clearly different (both in terms of carbon and nitrogen) across these same environmental contexts. Again, a good discussion here would be about why it might be beneficial for species to specialize based on the resources most optimal for them in a given environment.
CD <- lm(CD ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(CD)
## Fixed term is "(Intercept)"
Again, high \(\Delta\)AIC.
CD_best <- lm(CD ~ 1,
data = island_metrics)
sim_CD <- simulateResiduals(fittedModel = CD_best, plot = T)
This interpretation is really similar to that of TA, since these values are variations of the same measure of trophic space. Predators in different environmental contexts do not have broader trophic niches, even though their trophic niche is clearly different (both in terms of carbon and nitrogen) across these same environmental contexts. Again, a good discussion here would be about why it might be beneficial for species to specialize based on the resources most optimal for them in a given environment.
NND <- lm(NND ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(NND)
## Fixed term is "(Intercept)"
This \(\Delta\)AIC value is close, so we can look at that factor prod_level in a model for marginal significance.
NND_best <- lm(NND ~ prod_level,
data = island_metrics)
sim_nnd <- simulateResiduals(fittedModel = NND_best, plot = T)
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible
## Unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.
summary(NND_best)
##
## Call:
## lm(formula = NND ~ prod_level, data = island_metrics)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.02482 -0.13540 -0.01111 -0.04674 -0.02837 -0.01371 0.21050
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.50657 0.05769 8.781 0.000318 ***
## prod_levellow 0.17789 0.08812 2.019 0.099539 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1154 on 5 degrees of freedom
## Multiple R-squared: 0.449, Adjusted R-squared: 0.3388
## F-statistic: 4.075 on 1 and 5 DF, p-value: 0.09954
Yep, marginal significance.
plot(allEffects(NND_best))
There doesn’t seem to be a strong significant effect of environmental context on nearest neighbor distance. This means that predators tend to be trophically similar to each other across environmental contexts. We did see a marginal effect of island productivity on this value, with predators on low productivity islands having slightly higher values for nearest neighbor distance (\(\beta\) = 0.18, p-value = 0.10). What this means is that on low productivity islands, predators are slightly more trophically distinct within a population, maybe related to the need to search further for resources. Again, this is only a marginal effect.
SDNND <- lm(SDNND ~ prod_level*size_level,
data = island_metrics,
na.action = "na.fail")
dredge(SDNND)
## Fixed term is "(Intercept)"
Again, high \(\Delta\)AIC value.
SDNND_best <- lm(SDNND ~ 1,
data = island_metrics)
sim_sdnnd <- simulateResiduals(fittedModel = SDNND_best, plot = T)
There doesn’t seem to be a strong significant effect of environmental context on the standard deviation of nearest neighbor distance. This means that predators tend to be trophically similar to each other across environmental contexts. Given that this metric is less sensitive to sample size than NND - I would trust this result more than the marginal result above
size_is <- glmmTMB(Mass_g ~ prod_level*size_level + (1|Island),
data = spider_iso)
dredge(size_is)
## Fixed terms are "cond((Int))" and "disp((Int))"
\(\Delta\)AIC is within that 2AIC value cutoff, so I’ll go ahead and explore this likely marginal effect
size_is_best <- glmmTMB(Mass_g ~ size_level + (1|Island),
data = spider_iso)
size_is_sim <- simulateResiduals(fittedModel = size_is_best, plot=T)
summary(size_is_best)
## Family: gaussian ( identity )
## Formula: Mass_g ~ size_level + (1 | Island)
## Data: spider_iso
##
## AIC BIC logLik deviance df.resid
## 28.1 37.5 -10.1 20.1 72
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Island (Intercept) 0.02059 0.1435
## Residual 0.06666 0.2582
## Number of obs: 76, groups: Island, 7
##
## Dispersion estimate for gaussian family (sigma^2): 0.0667
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.36242 0.11609 3.122 0.0018 **
## size_levelsmall 0.06841 0.13720 0.499 0.6180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ran this anyway, but a pretty clearly non-significant difference based on island size in body size.
plot(allEffects(size_is_best))
#### Size by island interpretation
Even though predators are trophically distinct across high-low productivity islands (and also across big-small), this does not translate to individuals in the population reaching greater body sizes on islands with higher trophic positions. Although they eat more predators, predators on higher productivity islands aren’t larger. This could be evidence of the underlying increase in trophic complexity suggested in Young et al. 2013 on high productivity islands. A future exploration could be whether these predators, even though similar in size, could be different in population demographics (e.g population size, fecundity) across these environments based on resource availability.
C_size <- glmmTMB(d13C ~ Mass_g + (1|Island),
data = spider_iso)
dredge(C_size)
## Fixed terms are "cond((Int))" and "disp((Int))"
size_c_sim <- simulateResiduals(fittedModel = C_size, plot=T)
summary(C_size)
## Family: gaussian ( identity )
## Formula: d13C ~ Mass_g + (1 | Island)
## Data: spider_iso
##
## AIC BIC logLik deviance df.resid
## 190.2 199.5 -91.1 182.2 72
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Island (Intercept) 1.0253 1.0126
## Residual 0.4799 0.6928
## Number of obs: 76, groups: Island, 7
##
## Dispersion estimate for gaussian family (sigma^2): 0.48
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -24.3288 0.4125 -58.98 <2e-16 ***
## Mass_g 0.5989 0.3204 1.87 0.0616 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(C_size))
#### raw \(\delta^{13}C\) ~ size Interpretions
Larger individuals have a slightly higher \(\delta^{13}C\) value, suggesting that they are eating more marine, slightly.
N_size <- glmmTMB(d15N_c ~ Mass_g + (1|Island),
data = spider_iso)
dredge(N_size)
## Fixed terms are "cond((Int))" and "disp((Int))"
N_size_best <- glmmTMB(d15N_c ~ 1 + (1|Island),
data = spider_iso)
n_size_sim <- simulateResiduals(fittedModel = N_size_best, plot=T)
Body size does not seem to be determining trophic position measured by \(\delta^{15}N\). This could be because smaller predators can take advantage of smaller predators and that smaller species which can be consumed by smaller individuals may also be predators.
It may be a good idea to think about the SEAc values calculated by the groupMetricsML() function as a better value for determining uncertainty in measurements for capturing the values for both TA and CD described above in the metrics from Layman. Looks like I think this is the citation: https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.1365-2656.2011.01806.x
groupMetricsML() requires that you create a SIBER object like you did before, so maybe if you can find those values, you can also run a model on those similar to the linear models performed above.